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ABSTRACT 

Modern sky surveys are returning precision measurements of cosmological statistics such as weak lensing 
shear correlations, the distribution of galaxies, and cluster abundance. To fully exploit these observations, the- 
orists must provide predictions that are at least as accurate as the measurements, as well as robust estimates of 
systematic errors that are inherent to the modeling process . In the nonlinear regime of structure formation, this 
challenge can only be overcome by developing a large-scale, multi -physics simulation capability covering a range 
of cosmological models and astrophysical processes. As a first step to achieving this goal, we have recently de- 
veloped a prediction scheme for the matter power spectrum (a so-called emulator), accurate at the 1 % level out to 
k ~ 1 Mpc -1 and z = 1 for wCDM cosmologies based on a set of high-accuracy N-body simulations. It is highly 
desirable to increase the range in both redshift and wavenumber and to extend the reach in cosmological param- 
eter space. To make progress in this direction, while minimizing computational cost, we present a strategy that 
maximally re-uses the original simulations. We demonstrate improvement over the original spatial dynamic range 
by an order of magnitude, reaching k ~ 10 /zMpc" 1 , a four-fold increase in redshift coverage, to z = 4, and now 
include the Hubble parameter as a new independent variable. To further the range in k and z, a new set of nested 
simulations run at modest cost is added to the original set. The extension in h is performed by including pertur- 
bation theory results within a multi-scale procedure for building the emulator. This economical methodology still 
gives excellent error control, ~ 5 % near the edges of the domain of applicability of the emulator. A public domain 
code for the new emulator is released as part of the work presented in this paper. 

Subject headings: methods: statistical — cosmology: large-scale structure of the universe 



1. INTRODUCTION 
1.1. Background 

The discove r y of th e ac celerated expansi o n of t he Universe 
by Riess et alj (11998?) and iPerlmutter et al.l (Il999l) . now con- 
firmed by multiple observations, has had a tremendous impact 
on the field of observational and theoretical cosmology. Due to 
our near-complete ignorance regarding the origin of cosmic ac- 
celeration, the initial task is to better characterize the nature of 
the underlying cause, first by observational means. The primary 
focus in this endeavor is to decide whether the acceleration is 
caused by some form of dark energy, a cosmological constant 
being the simplest realization of this idea, or by some modifi- 
cation of general relativity on large cosmological scales. 

Unfortunately, carrying out controlled experiments in cos- 
mology is impossible - therefore, information relevant to the 
above task must be obtained by combining inferences from dif- 
ferent kinds of observations. These observations include the 
measurement of the temperature anisotropy of the cosmic mi- 
crowave background (CMB), the luminosity distance relation 
from supernova observations, and large scale structure probes, 
such as baryon acoustic oscillations (BAO), weak lensing, clus- 
tering of galaxies, and the abundance of galaxy clusters. De- 
spite the impressive power of the observations, all cosmological 
probes are plagued with a large number of sources of system- 
atic error. These can be broken up into two kinds, (i) those re- 
lated to the observational techniques - shape measurements in 
weak lensing and the determination of cluster masses being two 



prominent examples, and (ii) those related to the uncertainties 
in the theoretical predictions and the (related) errors in solving 
cosmological statistical inverse problems. An important mem- 
ber of this latter class is the prediction accuracy of the matter 
power spectrum, the central theme of this paper. 

Because robust results in cosmology must rely on combin- 
ing a number of observational probes, each possessing their 
individual strengths and weaknesses (not all of which can be 
fully predicted in advance), the current strategy is to proceed 
with several survey missions. In the spectroscopic realm, on- 
going or recently completed surveys focusing on the measure- 
ment of BAO and gala xy clustering inclu de the Sloan Digi- 
tal S ky Survey (SP SS) dZehavietalJ 120111) . the WiggleZ sur- 
vey dDrinkwater et a l. 2010) and the Baryon Oscillation Spec - 



troscopic Survey (BOSS) dSchlegel. White. & Eisensteinll2 009) 
component of SDSS-III; planned next-generation missions in- 
clude the Mid-Scale Dark Energy Spectroscopic Instrument 
(MS-DESI) and SuMIRe. Imaging surveys include SDSS, the 
Canada-France-Hawaii Telescope Legacy Survey (CFHTLS), 
the recently started Dark Energy Survey (DES), and in the fu- 
ture, Subaru Measurement of Images and Redshifts (SuMIRe) 
(Sugai 2012), the Large Synoptic Survey Telescope (LSST) 
(Ab elletaljl2009l) . and Euclid dRefregier et all 12010b . The 
most recent constraints from wea k lensing hav e been obtained 
by the CFH TLenS team (Kilbingeret al. 2012) and analysis of 
SDSS data dHuff et alJ lToi ll iLin et alj|201 1|). In the x -rav the 
next planned mission is eROSITA (M erloni et al. 2012), which 
will provide a significant update to results from prior surveys, 
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such as the 400d cluster survey (IVikhlinin et al.l 120091) . In the 
CMB, high-resolutio n measurements by the Atacama Cosmol- 
ogy telescope (ACT) (Das et al. 2012) and the South Pole Tele- 
scope (SPT) (van Engelen et al. 2012) have confirmed the sen- 
sitivity of the CMB to the large scale distribution of matter 
(CMB lensing). Significantly, firs t results from th e Planck mis- 
sion have been recently released (lAde et al.ll2013l) . For a com- 
prehensive review of the cu rrent sta t e-of-th e-art (pre-Planck) 
we refer the reader to Weinber g et al.l(t2012b . 

With the advent of ambitious large-scale observational pro- 
grams, emphasis on statistical errors will no longer be the sole 
driving force. The understanding, control, and reduction of 
systematic errors is the next major challenge - there are sev- 
eral cases where measurements are already systematics-limited. 
Additionally, in the nonlinear regime of structure formation, 
state of the art measurements are pushing up against the lim- 
its of available prediction accuracy from theoretical computa- 
tions and associated modeling. Major simulation campaigns, 
in concert with observational inputs, will be needed to control 
the error budgets. The simulation-based theoretical modeling 
process - consisting of a combination of first-principles calcu- 
lations and parameterized models - will have to account ac- 
curately, and robustly, for the nonlinear evolution of structure 
formation, gas physics, and complex feedback effects. 

Two examples will suffice to provide illustrations of the sort 
of difficulties that have to be faced. The first is cluster cos- 
molo gy - it was recently shown by ICunha & Evrardl ([2010) 
and IWu. Z entner. & Wechslei] (1201 Oh that the halo mass func- 
tion has to be predicted at the 1% level of accuracy for a DES- 
like survey in order to avoid errors on dark energy parameters. 
At the same time, an error in the prediction of halo masses 
at the 2% level will translate into mass functio n uncertainties 
in the cluster mass range of 5-10% (Cf. Bhattach arva et al. I 
2011). The effect is particularly strong because of the ex- 
ponential fall-off of the mass function at high masses. The 
second example relates to prediction error requirements for 
weak lensing shear. These are also severe - for the matter 
power spectrum, predictions at the percent level accuracy at 
k ~ 10 /r'Mpc are needed to extract all available information 
on dark energy (Huterer & Takada 2005). A more recent study 
by iHearin & Zentnerl (|2011|) concludes that the accuracy re- 
quirements are even more stringent - for imaging surveys such 
as LSST and Euclid, the power spectrum needs to be predicted 
at 0.5% accuracy out to k ~ 5 /i _1 Mpc. 

1.2. Fitting Functions vs. Emulators 

Currently, it is common practice to use fitting functions mo- 
tivated by theoretical heuristics and matched to simulations 
to provide the required modeling input. For instance, analy- 
ses of weak lensing measurements employ fitting functions for 
the matter power spectrum. ISem boloni et al. (2006) us ed the 
Peacock and Dodd s fitting function ([Peacock & Dodds 1996) 
as well as Halofit (Sm ith et al.ll2003l) f or their analysis of the 
CFHTLS cosmic shear measurements. iLin et al] (1201 1 ) used 
Halofi t for their recent SDSS analysis, as do Kilbi nger et al.l 
(2012) in the most recent CFHTLens analysis (they point out 
that an emulator would be more accurate but one was not avail- 
able over the required redshift range). Over a limited k range, 
the fitting functions are accurate at the ~ 5 — 10% level for 
ACDM models (see, e.g jHeitmann et al. 2010). The latest im - 
provement of Halofit is provided by Takahash i et al.l 120121: a 
comparison of our results with theirs is provided in Section [ 



For next-generation applications, the fitting function approach 
suffers from two deficiencies, first, error-control is non-uniform 
over cosmological model parameter space, and, second, the fit- 
ting process degrades the actual accuracy of the simulations 
used to build the parameterized fit. Moreover, because the fit- 
ting forms become essentially arbitrary as the accuracy require- 
ments become more stringent, this strategy is difficult - if not 
impossible - to implement systematically across a large number 
of freely floating cosmological and modeling parameters. 

Cluster mass functions provide another example of these dif- 
ficulties - obtaining dark energy constraints from the abun- 
dance of clusters of galaxies as carried out in, e.g., Vikhlinin et al. 
(2009), employs fitting forms for the mass function; to extend 
these fits across cosmological models, the assumption of uni- 
versality (Jenkins et al. 2001) is required. However, as shown 
in, e.g., iBhattacharva et al7l (1201 lb . universality for wCDM 
models only holds at the ~ 10% level of accuracy. It appears 
quite unlikely that a simple mass function fit - valid for a large 
class of dark energy models, covering a wide redshift range, 
and satisfying percent level accuracy requirements - will ever 
be attainable. 

We have recently developed the "Cosmic Calibration Frame- 
work" (CCF) to provid e accurate prediction schemes for cos- 
mological observables (Heit mann et al.l2006llHabib et al.l2 007) 
seeking to avoid the shortcomings of the fitting function ap- 
proach mentioned above. The CCF contains error-controlled 
direct numerical oracles for the predicted quantities, as opposed 
to (potentially uncontrolled) analytic fitting forms. It aims to 
make tools available that act as "LSSFast", essentially instan- 
taneous predictions of large scale structure observables, e.g., 
the nonlinear power spectrum, mass functions for different halo 
definitions, or the halo concentration-mass relation. 

At the heart of the CCF lies a sophisticated sampling scheme 
for optimally placing a given finite number of cosmological 
models in parameter space. Simulations are carried out at these 
parameter values (we use orthogonal-array Latin hybercubes as 
well as symmetric Latin hybercube designs). The next step is an 
efficient representation that translates the measurements from 
the simulations into functions that are conveniently interpolated 
(we use a principal component basis to provide a reduced data 
representation), and finally an accurate interpolation scheme 
over the basis functions (our choice here is Gaussian process 
mod eling). For an intro duction to the general framework, see, 
e.g. JSantner et al.ll2003l 

The CCF itself was first describ ed in Hei tmann et al. (2006), 
with details and examples given in lHabib et al] d2007l) . In a se- 
ries of three papers (Coyote Universe I-III), the development 
of a precision emulator for the matter power spectrum over a 
five-dimensional parameter space (8 = {u)b,w m ,n s ,w, <?%}) was 
described. This emulator provides predictions for the power 
spectrum for wCDM cosmologies out to k ~ 1 Mpc" 1 at the 
1% accuracy level for a redshift range of < z < 1, over 
a relatively wide range of cosmological parameters (see Sec- 
ti on [3]). Since then, similar appro ache s have been followed 
in lSchneider, Holm. & Knoxl d201 ll) and lAgarwal et al] J2012I) 
with some differences in the sampling and interpolation schemes 
used to build power spectrum emulators; 

The CCF framework was extended in lSchneider et al.l (12008) 
to derive an approximate statistical model for the sample vari- 
ance d istribution of the nonlinear matter power spectrum. Eifler 
(201 1 ) used the emulator to generate a weak-lensing prediction 
code to calculate various second-order cosmic shear statistics, 
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e.g., shear power spectrum, shear-shear correlation function, 
ring statistics and Complete Orthogonal Set of EB-mode In- 
tegrals (COSEBIs). The emulator presented in lLawrence et al.l 
d2010h was used bv lHuff etal1d2011l) in their analysis of SDSS 
weak lensing (where it was combined with Halofit to obtain the 
power spectrum at larger k values). More CCF-based emulators 
are under development; a recent e xample is an emul ator for the 
halo concentration-mass relation (Kwan et al. 2012D. 

Observational require ments, s uch a s tho se for PES weak 



lensing, and the work by lEiflerl (1201 ll) and iHuff et al.1 (1201 ll) 
have prompted us to extend the lLawrence et al.l (1201 Ol) power 



spectrum emulator in three directions: (i) making the Hubble 
parameter h freely choosable (in the original emulator the value 
for h in each model is fixed by the CMB distance to last scatter- 
ing constraint); (ii) extending the k range out to k = 8.6 Mpc -1 ; 
(iii) extending the z range out to z = 4. As described in Section[3] 
below, these extensions are nontrivial - a brute force approach 
would require simulations of size L = 1.5-2 Gpc with at least 
Np = 10,000 3 parti cles to fulfill the requirements for 1% accu- 
racy, as derived in iHeitmann et al.l (|2010|) . The most stringent 
requirement is due to the fact that the power spectrum should 
not be measured further out than half of the particle Nyquist 
frequency, k^ y = tt/A p with A ; , being the initial particle sepa- 
ration L/Np, and at the same time the linear modes should be 
well resolved in the simulation volume. In addition, at high red- 
shift and large k, particle shot noise becomes a significant issue 
(see the discussion in Section|3]l. 

Running a large number of the brute force simulations de- 
scribed above is currently impractical. In order to extend our 
prediction scheme we therefore choose another route, by us- 
ing a set of nested simulations of various volumes, and match- 
ing the results together at different values of k and z. The 
box lengths are chosen to be L = 1300 Mpc, L = 365 Mpc, 
L = 180 Mpc, and for the very high redshift results, L = 90 Mpc. 
This approach has its deficiencies, as discussed further in Sec- 
tion [3] Nevertheless, we can demonstrate that our predictions 
hold at the level of better than 5% error over the full k range. 

1 .3. The Power Spectrum at Small Scales: Baryonic Effects 

At small length scales, complex baryonic processes become 
important. Because these are difficult to model, a relatively 
large uncertainty in the power spectrum exists over the upper k 
range treated in t his paper. (Massive neutrinos a ls o have a small 
effect, see, e.g.. iBird. Viel. & Haehnelll l2012h IWhitel ri2004) 
estimated the effect of baryonic cooling on the mass power 
spectrum using the halo model, finding an increase in power 
at k ~ 10 /zMpc -1 a t the level of a few per cent at z = 0. An anal- 
ogous approach by Zhan & Knox ( |2004|) to study the effect of 
the intracluster medium (hot baryons) found a suppression in 
the power spectrum at roughly similar levels. 

Shortly after these papers appeared, simul ations were car- 
ried out to study these effects in more detail. IJing et al. (2006) 
carried out two simulations with the smoothed-particle hydro- 
dynamics (SPH) code Gadget-2, a non-radiative gas simula- 
tion and another including radiative cooling and star formation, 
with a third gravity-only simulation serving as a reference. At 
k ~ 10 hMpc~ l they found an effect on the total matter power 
spectrum at the few percent level in the non-radiative gas simu- 
lation and at the 10% level in the gas simulation with radiative 
cooling and star formation at z = (the effects are smaller at 
higher redshifts - for the radiative cooling/star formation sim- 
ulation, they are at the 2-3% level at z = 1). In both cases, the 



baryonic effects led to an enhancement in t he power spectrum . 
A similar campaign was later carried out bv lRudd et al.l (12008b : 
gravity-only, non-radiative gas dynamics, and a third simulation 
including radiative heating and cooling of baryons, star forma- 
tion, and supernova feedback. This set of simulations used the 
Adaptive Refinement Tree (ART) code, combining an N-body 
method and an Eulerian hydrodynamics solver on an adaptive 
mesh. Fo r the non-r adiative gas simulation they found a similar 
effect as IJing et all d200 6) on the overall matter power spec- 
trum, though effects on the baryon and dark matter components 
separately were different. The results from their third simula- 
tion showed a dramatic incr ease of the matte r power spectrum 
well beyond the effect that IJing et al.l (2006) observed (at the 
50% level at k ~ 5 /iMpc" 1 ) - but as noted by the authors, this 
result was not meant to be definitive. 

More recent studies have been carried out by Ivan Daalen et all 
(2011), again using Gadget-2 (the authors also provide a 
much more comprehensive discussion of results from other 
groups than we have space to do here). The main difference in 
this study is the inclusion of feedback from active galactic nu- 
clei (AGN). They observe a strong effect on the power spectrum 
in the opposite direction than found by the previous studies: a 
decrease in the matter power spectrum at 1 % at k ~ 0.3 /zMpc -1 , 
at 10% at k ~ 1 hMpc~\ and at 30% at k ~ 30 hMpc~ l . As 
shown by these results, obtaining accurate first principles pre- 
dictions for the power spectrum including baryonic effects will 
be very difficult. Multiple issues, ranging from physical effects 
to simulation uncertainties, remain to be sorted out. 

A strategy to properly deal with these difficulties is still 
emerging, but it will definitely require methods to incorporate 
information from observations, baryonic simulations, and an 
accurate calibration of the gravity-only power spectrum (which 
has no free modeling parameters). The difficulty of accounting 
for all "gastrophysics" effects at the desired accuracies, almost 
certainly precludes using hydro simulations in a fully predic- 
tive mode. One possible pathway is to construct parameterized 
models for baryonic effe cts on top of predictive N-body sim- 
ulatio n results (see e.g., Sembo loni et al.ll2011t IZentner et al.l 
2012 for some recent work) and combine these with (possibly 
multiple) observations ('self-calibration'), to successfully ex- 
tract cosmological information. In this paper, our target is to 
establish the bedrock on which all such analyses must rest - an 
accurate calibration of the gravity-only matter power spectrum. 

The rest of the paper is organized as follows. After provid- 
ing a brief summary of our power spectrum estimation from the 
simulations (Section |2), we describe the cosmological model 
space covered and the simulation suite used to build the emula- 
tor (Section|3]l. In Section|4]we outline our strategy to include 
the Hubble parameter h as a free parameter, without adding 
more N-body simulations. Next we discuss the generation of 
smoothed predictions for the power spectra for each model; this 
process underlies the interpolation scheme for constructing the 
emulator. In Section|6]we show some examples from the work- 
ing emulator, including a brief comparison with Halofit. Fi- 
nally, we end with a conclusion and outlook in Section [7] The 
Appendix n presents a short discussion of an improvement to the 
general Gaussian process approach employed here compared to 
that in the previous Coyote papers, including a minor numerical 
correction of the earlier results. 

2. POWER SPECTRUM ESTIMATION 

The key statistical observable in this paper is the density fluc- 
tuation power spectrum P(k), the Fourier transform of the two- 
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point density correlation function. We follow the same strat- 
egy for extracting the power spectra from the simulations as in 
Heitman n et alj d2010l) and give a summary here for complete- 



ness. 



In dimensionless form, the power spectrum may be written 



as 



A 2 (k) = 



k 3 P(k) 



(1) 



which is the contribution to the variance of the density pertur- 
bations per lnfc. 

Because A^-body simulations use particles, one does not di- 
rectly compute P(k) or equivalently, A 2 (k). Our procedure is to 
first define a density field on a grid of size iV? with a fine enough 
resolution such that the grid filtering scale is much higher than 
the k scale of interest. This particle deposition step is carried 
out using Cloud-in-Cell (CIC) assignment. The application of 
a discrete Fourier transform (FFT) then yields <5(k) from which 
we can compute P(k) = <5(k)| 2 , which in turn can be binned 
in amplitudes to finally obtain P(k). Since the CIC assignment 
scheme is in effect a spatial filter, the smoothing can be com- 
pensated by dividing P(k) by W 2 (k), where 

'k x N g \ . 2 {kyN g \ .zfkyNg" 



W(k)=f 



(2) 



Typically the effect of this correction is only felt close to the 
maximum (Nyquist) wavenumber for the corresponding choice 
of grid size. One should also keep in mind that particle noise 
and aliasing artifacts can arise due to the finite number of par- 
ticles used in A^-body simulations and due to the finite grid size 
which is used for the power spectrum estimation, as discussed 
further in Section 13.11 For an extensi ve suite of convergen ce 
tests addressing these issues, see, e.g. jHeitmann et alj(l2010h . 
We average P(k) in bins linearly spaced in k of width Ak ~ 
0.001 Mpc -1 , and report this average for each bin containing at 
least one grid point. We assign to each bin the k associated with 
the unweighted average of the k's for each grid point in the bin. 
Note that this procedure introduces a bias in principle, since 
for nonlinear functions (f(x)) / f((x)), but our bins are small 
enough to render this bias negligible. Throughout the paper we 
will show results for both A 2 (k) and P(k). The emulator itself 
provides results for both definitions as well. 

3. COSMOLOGICAL MODELS AND SIMULATION SETS 

The emulator is based on 37 cosmological models spanning 
the class of wCDM cosmologies. We allow for variations of the 
following six parameters: 

9 = {uj h ,u> m ,n s ,h,w,a & }. (3) 

The 37 models are chosen to lie within the ranges: 

0.0215 <w/,< 0.0235, 

0.120 <lu,„< 0.155, 

0.85 <n s < 1.05, 

0.55 <h< 0.85, 

-1.30<w<-0.70, 

0.616 <ct 8 <0.9, 

motivate d by recent constrai nts from CMB measurements by 
WMAP JKomatsu et alj|201ll) . In our original work we locked 
the value of the Hubble parameter h to the best-fit value given 
a measurement of the distance to the s urface of large scattering 
for each model (lLawrence et alJ l2010). The values for h in the 
original runs then ranged from 0.55 < h < 0.85. While this 
range is larger than the currently available constraints on the 



(4) 



Hubble parameter from, e.g., iRiess et all (J2011) we decide to 
keep it in our new runs in order to include the best fit models. 
In Section H] we explain how we extend the parameter range to 
include h as a free variable in the new emulator without actually 
running more N-body simulations. In addition to the 37 models, 
we ran one ACDM model (M000 in Table[T]) which is not used 
to build the emulator, but is instead used as a reference to test 
the emulator accuracy. All 37+1 models are specified in detail 
in Tabled] 

The specific model selection proc ess for the Coyo t e Uni - 
verse runs is described at length in iHeitmann et al.l (|2009). 
It is based on Symmetric Latin Hypercube (SLH) sampling 
dLi & Yd l2000). Following the SLH strategy guarantees good 
coverage of the parameter hypercube. In our specific case we 
chose an SLH design that has space filling properties in the case 
of two-dimensional projections in parameter space. In other 
words, if any two parameters are shown in a plane, the plane 
will be well covered by simulation points. A n extensive dis- 
cussion of optimal design choices is given in Heit mann et al.l 
(2009). The interested reader is referred to that paper for de- 
tails. 

The emulator developed here is valid over the redshift range 
z= {0,4} and the Grange extends out to k= 8.6 Mpc" 1 . In units 
of /iMpc" 1 this covers k ranges between k ~ 10 /iMpc" 1 and k ~ 
15 /iMpc" 1 , depending on the particular cosmological model. 
In order to enable a smooth interpolation between redshifts, we 
store results at 1 1 outputs: 

a= {1.0,0.9, 0.8, 0.7, 0.6, 0.5,0.4, 

0.3333, 0.2857, 0.25, 0.2}. (5) 

We use different box sizes to cover different ranges in k and 
redshift space. A summary of the simulation sizes is given in 
Table|2] The largest set of simulat ions is from the original Coy- 
ote Universe suite as described in lLawrence et al.l (120101) . This 
set of runs (all carried out in a 1300 Mpc box) consists of - for 
each cosmological model - 16 realizations with 512 3 particles 
using a particle mesh (PM) code with a 1024 3 grid, 4 PM real- 
izations with 1024 3 particles run on a 2048 3 grid, a nd one high 
resolution GADGET-2 TreePM run (Springel 20051) with 1024 3 
particles. In addition, for the results reported here we run, for 
each model, simulations with 512 3 particles and 10 kpc force 
resolution with Gadget-2 in a 365 Mpc box (one realization 
per model), a 180 Mpc box (three realizations per model), and 
a 90 Mpc box (4 realizations per model). As we explain below, 
we do not run the 90 Mpc box to z = but stop at z = 2. We will 
discuss the reasons for our specific choices and how we match 
the different boxes in the following section. 

3.1. Nested Simulations 

The extension of the original emulator beyond k = 1 Mpc -1 
and z = 1 at high accuracy is nontrivial due to two limiting fac- 
tors: the spatial Nyquist frequency, 

; nN P «* 

fc Ny=— — , (6) 

setting the largest viable £-value, and the particle shot noise 
limit: 

P^{k)=(jf\ , (7) 

restricting the lowest amplitudes at which the power spectrum 
can be measured accurately. These iss ues have been k nown for 
a long time (see, e.g. lBaugh. Gaztanaga. & E fstathiou 1995] for 
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Table 1 
Parameters for the 37+1 models which define the sample space; k„i is measured in Mpc~'. See text for further details. 



# 


w m 


u> b 


n s 


— w 


CT8 


h 


# 


Wm 


UJ b 


n s 


— w 


cs 


h 


M000 


0.1296 


0.0224 


0.9700 


1.000 


0.8000 


0.7200 


M019 


0.1279 


0.0232 


0.8629 


1.184 


0.6159 


0.8120 


M001 


0.1539 


0.0231 


0.9468 


0.816 


0.8161 


0.5977 


M020 


0.1290 


0.0220 


1.0242 


0.797 


0.7972 


0.6442 


M002 


0.1460 


0.0227 


0.8952 


0.758 


0.8548 


0.5970 


M021 


0.1335 


0.0221 


1.0371 


1.165 


0.6563 


0.7601 


M003 


0.1324 


0.0235 


0.9984 


0.874 


0.8484 


0.6763 


M022 


0.1505 


0.0225 


1.0500 


1.107 


0.7678 


0.6736 


M004 


0.1381 


0.0227 


0.9339 


1.087 


0.7000 


0.7204 


M023 


0.1211 


0.0220 


0.9016 


1.261 


0.6664 


0.8694 


M005 


0.1358 


0.0216 


0.9726 


1.242 


0.8226 


0.7669 


M024 


0.1302 


0.0226 


0.9532 


1.300 


0.6644 


0.8380 


M006 


0.1516 


0.0229 


0.9145 


1.223 


0.6705 


0.7040 


M025 


0.1494 


0.0217 


1.0113 


0.719 


0.7398 


0.5724 


M007 


0.1268 


0.0223 


0.9210 


0.700 


0.7474 


0.6189 


M026 


0.1347 


0.0232 


0.9081 


0.952 


0.7995 


0.6931 


M008 


0.1448 


0.0223 


0.9855 


1.203 


0.8090 


0.7218 


M027 


0.1369 


0.0224 


0.8500 


0.836 


0.7111 


0.6387 


M009 


0.1392 


0.0234 


0.9790 


0.739 


0.6692 


0.6127 


M028 


0.1527 


0.0222 


0.8694 


0.932 


0.8068 


0.6189 


M010 


0.1403 


0.0218 


0.8565 


0.990 


0.7556 


0.6695 


M029 


0.1256 


0.0228 


1.0435 


0.913 


0.7087 


0.7067 


M011 


0.1437 


0.0234 


0.8823 


1.126 


0.7276 


0.7177 


M030 


0.1234 


0.0230 


0.8758 


0.777 


0.6739 


0.6626 


M012 


0.1223 


0.0225 


1.0048 


0.971 


0.6271 


0.7396 


M031 


0.1550 


0.0219 


0.9919 


1.068 


0.7041 


0.6394 


MOB 


0.1482 


0.0221 


0.9597 


0.855 


0.6508 


0.6107 


M032 


0.1200 


0.0229 


0.9661 


1.048 


0.7556 


0.7901 


M014 


0.1471 


0.0233 


1.0306 


1.010 


0.7075 


0.6688 


M033 


0.1399 


0.0225 


1.0407 


1.147 


0.8645 


0.7286 


M015 


0.1415 


0.0230 


1.0177 


1.281 


0.7692 


0.7737 


M034 


0.1497 


0.0227 


0.9239 


1.000 


0.8734 


0.6510 


M016 


0.1245 


0.0218 


0.9403 


1.145 


0.7437 


0.7929 


M035 


0.1485 


0.0221 


0.9604 


0.853 


0.8822 


0.6100 


M017 


0.1426 


0.0215 


0.9274 


0.893 


0.6865 


0.6305 


M036 


0.1216 


0.0233 


0.9387 


0.706 


0.8911 


0.6421 


M018 


0.1313 


0.0216 


0.8887 


1.029 


0.6440 


0.7136 


M037 


0.1495 


0.0228 


1.0233 


1.294 


0.9000 


0.7313 



Table 2 
Box sizes, particle numbers, force resolution, and corresponding values for shot noise limit, Nyquist frequency, 

AND MASS RESOLUTION. 



Length 


A/ 3 
p 


Number of 


Force res. 


Shot noise 


^Ny 


m p 


Scale factors 


[Mpc] 




Realizations 


[kpc] 


[Mpc 3 ] 


[Mpc" 1 ] 


[M ] 




1300 


512 3 


16 


1270 


2.05 


2.48 


5.7 -10 n u; m 


0.2-1.0 


1300 


1024 3 


3 


635 


2.05 


2.48 


5.7 -lO'Vn 


0.2-1.0 


1300 


1024 3 


1 


50 


2.05 


2.48 


5.7 • 10 n u; m 


0.2-1.0 


365 


512 3 


2 


10 


0.36 


4.41 


1.0-10 n w m 


0.4-1.0 


180 


512 3 


3 


10 


0.04 


8.94 


1.2-10 10 w ffl 


0.2-0.6 


90 


512 3 


4 


10 


0.005 


17.87 


1.5-10 9 w m 


0.2-0.33 
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an early discussion). With the recent efforts to obtain very ac- 
curate results for the absolute measurement of the power spec- 
trum out to small scales, it is worthwhile to briefly revisit the 
essential arguments. 

Both the Nyquist limit and the shot noise level depend on the 
mass resolution and the size of the simulation volume and in 
both cases reasonable results can only be obtained for a large 
number of particles. Computational limits imply a necessarily 
finite number of particles; the obvious option is to shrink the 
simulation volume. But this option has its own pitfalls - lack of 
long-range power and increased sampling variance, for exam- 
ple - and will break down at some point as discussed quantita- 
tively below. 

The shot noise effect is particularly annoying at early times; 
at late times, once the power spectrum has risen substantially 
above the shot noise level, it becomes much less problematic. 
The Nyquist sampling issue, which essentially sets the dynamic 
range at which the initial condition can be produced, is worse 
for larger boxes with the same particle number. In addition, the 
impact of both problems depends on erg. For two simulations 
that differ only in their value of erg, the one with the higher erg 
will contain more nonlinear structure and growth at the same 
redshift - this means that the shot noise problem is worse for 
low erg runs while the issues with small volumes at lower red- 
shift are more severe for high erg models. 

At the current levels of available computational power, it is 
impossible to carry out a brute force approach to obtain a sub- 
percent accurate answer for the po wer spectrum at the de sired 
^-values. As discussed at length in Heitma nn et alj yOlO), the 
maximum value is set by k^ y /2. At the same time, in order to 
ensure that the largest modes in the box at z = are linear to 
high accuracy, a linear box size of 2000 Mpc would be a good 
choice. To reach a wavenumber value of 10 Mpc" 1 would imply 
a particle loading of Nt ~ 12700 3 ~ 2 trillion part icles. While 
it is p ossible to carry out a few such simulations (lHabib et al.l 
|2012|) . a large suite of them is clearly currently out of reach. 

This notional set-up would lead to a shot noise level of 
~ 4 • 10~ 4 . The model with the lowest erg in our simulation 
set is M019 with erg = 0.6159. The amplitude of the power 
spectrum at z = 4 is P(k = 10, z = 4) = 0.025 Mpc 3 , translating 
to ~ 60 times above the shot noise level and is, hence, safe. 
If we wanted to decrease the number of particles by choosing 
a smaller volume, 1000 Mpc would be barely large enough. 
In this case, ~ 250 billion particles would be needed to push 
out the Nyquist limit far enough. This would lead to a shot 
noise level of ~ 4 • 10~ 3 , dangerously only a factor of 6.25 be- 
low the amplitude of the power spectrum at k = 10 Mpc" 1 and 
Z = 4. Moreover, even this simulation size constitutes a cur- 
rently prohibitive expense if a full suite of simulations has to be 
performed. 

Due to these obstacles, a strategy for mixing and matching 
simulation sizes - differently for different redshifts - needs to 
be employed. Smaller simulation volumes are required for high 
redshift results while larger volumes will be used for lower red- 
shifts. Although such a procedure is effective as we show be- 
low, because of the uncertainties induced by having to sew to- 
gether multiple simulation results, and because the number of 
cosmological models used is still the same as in the original 
set, sub-percent accurate power spectra are not obtained over 
the new, much wider dynamic range in wavenumber and red- 
shift. At the small spatial scales that the current emulator goes 
to, however, current uncertainties in characterizing bary onic ef- 




FlG. 1 . — Comparison of the simulations with resummed perturbation theory 
for M037 (averaged over 20 realizations) at three redshifts. The black hori- 
zontal lines show the 1% error limits. The red lines show the matching points 
we choose for connecting perturbation theory to the simulation results. For 
1.0 < a < 0.5 the matching value is at k = 0.03 Mpc -1 and for 0.4 < a < 0.2 
it is at k = 0.1 Mpc -1 . 



fects clearly outweigh its inaccuracies. Therefore, over the ex- 
tended dynamic range, the constraints on the accuracy of the 
dark matter power spectrum are not uniform, and not as strin- 
gent near the limits of the range probed in wavenumber and in 
redshift. 

We quantify below some of the errors due to the fact that the 
power spectra at all k and z do not arise from a single overarch- 
ing simulation. Over its full range, and across all cosmological 
models, the emulator is nevertheless accurate to better than 5%. 
We now describe the details of our matching strategy for the 
power spectra at each redshift. 

For the largest scales, we u se resummed perturbation theory 
(RPT) using the Copter code (ICarlson. White. & Padmanabhanl 
2009); the code has been modified to allow for w ^-1 . We start 
by breaking up the models into three groups, depending on the 
power spectrum normalization, erg: 

erg < 0.7, Assembly 1, (8) 

0.7 < erg < 0.8, Assembly 2, (9) 

erg > 0.8, Assembly 3. (10) 

We use RPT for all three cases in the same way up to the fol- 
lowing k values: 

1.0<a<0.5: k = 0.03 Mpc" 1 

0.4<a<0.2: £ = 0.1 Mpc" 1 . 

The first cut is the same as was used in Lawrenc e et al.l (12010). 
while the second is more aggressive. Since the power spec- 
trum is more linear at higher redshifts, perturbation theory will 
be valid out to higher k. We exhibit the accuracy of RPT at 
the matching scales for the model with the highest erg, M037 
(erg = 0.9) - the most difficult case - in Fig. Q] The ratio of 
RPT is shown with respect to the simulation output (the aver- 
age from our 20 realizations with L = 1300 Mpc) at three dif- 
ferent redshifts, z = 0, 1.5,4. The red vertical line shows our 
matching point for perturbation theory. In all cases, the er- 
ror is below 1%, in very good agreement with the findings of 
ICarlson. White. & PadmanabhanT(l2009l) . 

Next, we must determine the maximum k values out to which 
to use the results from each of the different volumes. The indi- 
vidual matching strategies for the three assemblies depending 
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Table 3 
Assembly 1, <r 8 > 0.8 

_a RPT 1300 Mpc 365 Mpc 180 Mpc 90Mpc 

0.2 k<0.l 0.1 < A; < 0.4 0.4 < k < 1.0 1 < fe < 8.6 

0.25 k<0.l 0.KK0.4 0.4 < k < 1.0 1<A:<8.6 

0.2857 k<0.l 0.KK0.4 0.4 < k < 1.0 1 < fe < 8.6 

0.3333 k<0.l 0.1 < A; < 0.4 0.4<fe<1.0 l<fc<8.6 

0.4 k<0.1 0.03 <£< 0.4 0.4 < A: < 8.6 

0.5 A < 0.03 0.03 <k< 1.0 1.0< A< 8.6 

0.6 A < 0.03 0.03 < A < 1.0 1.0 < A; < 8.6 

0.7 A < 0.03 0.03 < k < 1.0 1.0 < A < 8.6 

0.8 A < 0.03 0.03 < k < 1.0 1.0 < A: < 8.6 

0.9 A < 0.03 0.03 <k< 1.0 1.0 < A < 8.6 

1.0 A < 0.03 0.03 < k < 1.0 1.0 < A: < 8.6 

Table 4 
Assembly 2, 0.7 < a s < 0.8 

_a RPT 1300 Mpc 365 Mpc 180 Mpc 90 Mpc 

"02 k<0.1 0.1 <k<0A 0.4<k< 1.0 1 < A: < 8.6 

0.25 k<0.l 0.KK0.4 0.4 < A: < 1.0 1<A<8.6 

0.2857 k<0.l 0.1 < A; < 0.4 0.4 < k < 1.0 1 < k < 8.6 

0.3333 fe<0.1 0.1 < A; < 0.4 0.4 < k < 1.0 1 < k < 8.6 

0.4 k<0.l 0.03 < A < 0.4 0.4 < A: < 8.6 

0.5 A < 0.03 0.03 <k< 1.0 1.0 < fe< 8.6 

0.6 A < 0.03 0.03 <A< 1.0 1.0 < A; < 8.6 

0.7 A < 0.03 0.03 < it < 1.0 1.0< A:< 8.6 

0.8 A < 0.03 0.03 < k < 1.0 1.0 < A; < 8.6 

0.9 A < 0.03 0.03 < k < 1.0 1.0 < it < 8.6 

1.0 k<0.03 0.03 < k < 1.0 1.0 < k < 8.6 



Table 5 
Assembly 3, o-g < 0.7 

_a RPT 1300 Mpc 365 Mpc 180 Mpc 90 Mpc 

"02 A: < 0.1 0.1 <A<0.4 0.4<A< 1.0 1 < A: < 8.6 

0.25 A<0.1 0.KK0.4 0.4 < k < 1.0 1 < A < 8.6 

0.2857 A<0.1 0.KK0.4 0.4 < A: < 1.0 1<A<8.6 

0.3333 A<0.1 0.1<A<0.4 0.4 < k < 1.0 1 < A < 8.6 

0.4 A<0.1 0.03<A<0.4 0.4<A<8.6 

0.5 A < 0.03 0.03 <A< 1.0 1.0<A<8.6 

0.6 A < 0.03 0.03 < k < 1.0 1.0<A<8.6 

0.7 A < 0.03 0.03 <k< 1.0 1.0<A<8.6 

0.8 A < 0.03 0.03 <A< 1.0 1.0 < A; < 8.6 

0.9 A < 0.03 0.03 < k < 1.0 1.0<A<8.6 

1.0 A < 0.03 0.03 < k < 1.0 1.0 < k < 8.6 
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on erg are summarized in Tables [3]-[5] As mentioned above, the 
matching point depends on the Nyquist wavenumber and, at 
higher redshifts, the shot noise level. As previously established 
in Heitmann et al. (2010) the results from the 1300 Mpc volume 
hold at the sub-percent level accuracy out to k ~ 1 Mpc" 1 be- 
tween < z < 1 or 0.5 < a < 1.0. Therefore, the power spectra 
for all three assemblies in this range are taken from these sim- 
ulations (the simulations are also particularly valuable because 
of the large numbers of realizations for each of the boxes). 

At higher redshifts, the shot noise level of the earlier simu- 
lations becomes too high and we can use them only to a lower 
lvalue, consequently, in all three assemblies, for 0.2 < a < 0.4 
we cut off the power spectra from the 1300 Mpc boxes at 
k = 0.4 Mpc" 1 . For a > 0.5, the upper k value is taken to be 
1.0 Mpc" 1 . The following step is to match the results from the 
smaller boxes at these cutoff points. Shot noise level consider- 
ations play the dominant role in choosing the matching points 
for the power spectra from the different simulations. 

In the case of Assembly 1 (high erg), the second-largest 
box (365 Mpc) is used to cover the remaining full range in 
wavenumber from 0.4 < a < 1. The next two smaller boxes of 
size 180 Mpc and 90 Mpc are used to fill in the higher k range 
as shown in Table[3] The results of analogous strategies for As- 
sembly 2 (medium erg) and Assembly 3 (low erg) are given in 
Table EJand Table El 

Figure |2] shows the matched power spectra for two example 
models, M0 19 and M037, representing the lowest and high- 
est erg in our model selection. The upper panels show the full 
power spectra from all boxes while the lower panels show the 
results with cut-offs applied at the appropriate fc-values. The re- 
sults in the lower panel demonstrate that the approach of using 
nested volumes is quantitatively satisfactory. In the next two 
sub-sections we discuss the main errors involved in the match- 
ing procedure, (i) finite mass resolution which leads us to push 
beyond k^ y /2 for the small boxes, and (ii) finite volume effects. 

3.2. Finite mass resolution effects 

In iHeitmann et alj ([2010) the effects of finite mass resolu- 
tion on power spectrum estimation were discussed in some 
detail. Two effects investigated there act in opposite direc- 
tions. The first is sampling noise, i.e., the effect on computing 
a power spectrum from a density field computed from a dif- 
ferent number of particles. For this, a particle distribution at 
z = was down-sampled by factors of 8, 64, and 512 and the 
power spectra re-measured. The effect of the insufficient sam- 
pling of the density field led to an enhancement of the power 
spectrum larger than 1% beyond k^ y /2. The second is the ef- 
fect of lower particle sampling in the initial conditions. This 
leads to a suppression of the power spectrum since the initial 
conditions now lack small-scale power. This effect decreases 
for simulations with larger Nyquist frequency since the num- 
be r of modes in the box at small scales increases. In the test 
in IHeitmann et alj d2010h it was shown that a simulation with 
£Ny=0.859 /iMpc" 1 has lost more than 10% of power at k^ y /2, 
while for £ Ny =3.44 /iMpc" 1 it is well below 1% at k Ny /2. 

In the current paper, we use the large volume simulations 
(1300 Mpc) only up to k < 1 Mpc" 1 at low redshift and k < 
0.4 Mpc" 1 at high redshift, both values being well below k^ y /2. 
For the smallest volumes (90 Mpc) we also do not use any re- 
sults beyond k^ y /2, simply because for the smallest boxes the 
Nyquist limit is much larger than k ~ 10 Mpc" 1 . For the inter- 
mediate boxes (180 Mpc and 365 Mpc) we choose to be more 



aggressive. For the lowest redshifts we push the 365 Mpc limit 
to twice the Nyquist criterion and for intermediate redshifts we 
use the 180 Mpc box out to the Nyquist frequency. In princi- 
ple, we could have avoided pushing the 365 Mpc box beyond 
the Nyquist limit by using the 180 Mpc box at large k values 
instead. There are two reasons why we chose not to do this. 

First, the comparison of inaccuracies due to finite box effects 
versus going beyond the Nyquist limit suggests that the finite 
box errors dominate. The quasi-linear regime in the small box 
simulations close to the final time step at which they are used 
is not accurate because of the missing large-scale power. In ad- 
dition, the scatter in the overall amplitude is larger due to the 
small box size as discussed in the next sub-section. These two 
errors combine to outweigh the inaccuracy due to the Nyquist 
limit. This point is illustrated in Fig. [2] In the right upper 
panel (M037), at z = 0, the green curve shows the result for 
the 365 Mpc box while the blue curve shows the result for the 
180 Mpc box. It is apparent that the small box clearly does not 
capture the nonlinear turn-over behavior beyond k ~ 0.1 Mpc" 1 
very well. On the other hand, both results are in close agree- 
ment at large k, so the accuracy of the 365 Mpc box results is 
still good at that point. 

Second, the matching procedure described in Section |5] for 
melding the different power spectra pieces induces another 
small inaccuracy. This is again mainly due to finite box size 
effects, as the amplitude has to be adjusted to enable a match- 
ing of the boxes in the high k region. It is therefore desirable to 
minimize the number of matching points, and hence to take the 
results from the 365 Mpc boxes out to higher k values, rather 
than to introduce another matching point at intermediate scales. 
In order to verify that using the 365 Mpc box at higher k does 
not introduce a major error, we also checked that the difference 
at high k between the 365 Mpc box and the 180 Mpc box is 
small, as can be seen in Fig. [2] For M037 (high erg) we find 
differences of 2% and less beyond k = 6 Mpc" 1 and for M0 19 
(low erg) we find differences well below 1% in the high k range. 

To summarize this discussion, for the high redshift results, 
we do not use simulations beyond k^ y /2 (see Tables [3]|5] for 
more details) and therefore estimate the error throughout the 
k range as being well below 1%, based on the test results from 
Heit mann et alJJ2010h . For the low redshift cases, we use some 
results that go beyond the Nyquist limit. By comparing to the 
smaller boxes we estimate that the error does not exceed 2%. 

3.3. Finite box size effects 

As outlined in the Introduction, ultimately one would want 
to carry out simulations in large box sizes, 1-2 Gpc, but this is 
currently impractical for a large suite of runs, each with suf- 
ficient mass resolution. Smaller simulation volumes have two 
drawbacks: (i) At redshifts close to z = the undersampling 
of large-scale power may be a problem (the precise extent de- 
pending on erg and the box size). In smaller boxes, while the 
largest-scale modes might appear to be linear, their amplitude 
is actually suppressed by missing nonlinear power (in a large 
volume run, the same modes would have higher amplitudes), 
(ii) The realization scatter in small volumes is much larger 
due to the smaller number of large sc ale modes. The first o f 
these points was considered in depth in IHeitmann et alj (|2010). 
There, we demonstrated the suppression of the power spectrum 
due to finite volume effects by comparing the average of 4 real- 
izations in a 2000 /r'Mpc box, 8 realizations in a 936 /t _1 Mpc 
box, and 127 realizations in a a 234 /i _1 Mpc box. Figure 6 
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FIG. 2. — Two examples of how nested simulation volumes are used to cover a large range in k. Models shown have the most extreme values of erg, M019 
with <t 8 = 0.6159 (right column) and M037 with erg = 0.9 (left column). The upper row shows the full power spectra from the different simulations as well as the 
corresponding Nyquist limit (vertical line) and the shot noise limit (horizontal line). The lower row shows the different power spectra matched up (without any 
additional smoothing) at the /r-values described in Tables|3]-|5] The power spectra results are shown at z = and z = 4. In all cases the matching leads to a relatively 
smooth power spectrum covering the full k range of interest. In the lower panel, at z = 0, we show the good agreement between the 365 Mpc and 180 Mpc boxes at 
high k by overlaying the two power spectra (see text for details). Note that in the lower plots at z = we used the original emulator predictions for the intermediate 
k range shown in red. 



in lHeitmann et all d2010l) shows the suppression at mildly non- 
linear scales. The realization scatter was not a significant is- 
sue in that paper since the final simulation volumes were large 
(1300 Mpc) and we averaged over 20 realizations. 

In the current paper, the realization scatter is of concern since 
we do resort to using small volumes. The mildly nonlinear 
regime is still very accurate since it is covered by large volume 
simulations, so the key question is to understand the high-£ be- 
havior. In order to investigate this effect, we car ry out two tests. 
First, w e use the same 127 simulations used in lHeitmann et al.1 
(2010) to measure the overall dispersion in the power spectra 
between different realizations. These simulations were carried 
out in 325 Mpc volumes with a PM code almost matching the 
small simulation volumes used here. Since we are only inter- 
ested in the relative effect, these lower resolution simulations 
are sufficient. The results are encapsulated in Fig. [3] The upper 
panel shows the ratio of the average of the 127 power spectra 
with respect to each individual spectrum; in the lower panel, we 
attempt to correct the run-to-run scatter by simply adjusting the 
amplitude of each power spectrum to the average value at the 
highest measured k (which is well-sampled in each simulation). 
This procedure improves the scatter somewhat but is rather un- 
controlled and we decided not to use it. Nevertheless, the result 



is somewhat interesting. Overall, the run-to-run scatter at small 
scales is up to about 5%, at the matching scale of k ~ 1 Mpc" 1 , 
and up to 10% for the most extreme outliers. In order to reduce 
this problem, we carry out at least two realizations for each 
model, in some case up to four to obtain a realization in which 
the matching to the 1300 Mpc box is as close as possible. Since 
we match the even smaller volumes (180 Mpc and 90 Mpc) at 
larger k values, the problem there is not quite as severe. Nev- 
ertheless, even in those cases we run up to four realizations per 
model to provide an accurate match to the larger boxes. 

As a second check we carry out six simulations of the same 
cosmology each with box length L = 365 Mpc and A^ = 512 3 
particles, which are chosen to exactly match the size and resolu- 
tion of the smaller Coyote runs. The initial particle distribution 
is set with the same realization at z = 21 1 on large scales for all 
six simulations, but on scales smaller than k = 1.39 Mpc" 1 , we 
allow the initial conditions to vary between runs by changing 
the random seed for each simulation. These are then evolved 
using Gadget-2 and we show the ratio of power spectra at 
z = with reference to a single simulation in the set in Fig. [4] 
The red line marks where the initial power spectrum differs be- 
tween the realizations. Notice that the scatter in the power spec- 
trum is less than ^3% and appears unbiased. This gives an es- 
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FIG. 3. — Variations in power spectra due to different realizations from 
127 simulations, each for a 325 Mpc box size. The upper panel shows the 
ratio of each of the 127 results at z = with respect to the average power 
spectrum. The lower panel shows an attempt to correct the results - here we 
adjust the amplitude of each power spectrum so that it matches the average 
power spectrum at the largest i-value. While at large k this reduces the error 
somewhat, the procedure is not satisfactory overall and we do not use it in the 
final results. The scatter is roughly at the few percent level. 




rameters to be the best-fit CMB value. To allow for more flexi- 
bility, especially keeping in mind possible tensions in different 
observational values of h, we now aim to allow it to vary within 
the range that is covered by the original Coyote Universe runs, 
i.e., 0.55 < h < 0.85. This means that the sampled parameter 
space must be suitably extended. 

A first idea might be to use the existing 37 models (all of 
which have different values of h) and rebuild the emulator keep- 
ing the Hubble par ameter free. This nai ve approach is inade- 
quate. Figure 6 in lHeitmann et alj d2009l) (seventh row) shows 
the distribution of h values with respect to the other parame- 
ters. From that figure it is clear that the parameter space for h 
is not well covered if we restrict ourselves to the 37 available 
models. The design has large holes, particularly for high val- 
ues ofh. Our tests confirmed that this sub-optimal design does 
not lead to an accurate emulator. We note that this result is in 
fact a successful demonstration of our overall approach, since 
it shows that the optimality of the sampling design is indeed a 
crucial factor. 

The obvious alternate approach to include h as a new param- 
eter is to generate a new design over six parameters and in- 
crease the number of models to obtain sufficient accuracy. This 
is obviously undesirable for just one additional parameter since 
it would add a large computational cost and would not make 
use of the already available simulations. We therefore choose a 
quite different path. 

The new idea implemented here is to exploit a certain flex- 
ibility in constructing emulators; this flexibility relates to the 
fact that emulators can be constructed by incorporating results 
from multiple models across different scales, i.e., where the re- 
sult for each model does not necessarily have information avail- 
able across the complete set of scales. 

We proceed by generating 100 new predictions over a six- 
dimensional parameter space 9 = {LUt,,u) m ,n s ,h,w,cr$} withRPT 
out to k < 0.1 Mpc" 1 for high redshift (z > 2) and out to 
k < 0.03 Mpc" 1 for low redshift (z < 2). For larger fc-values we 
use the results from the original 37 simulations for which we al- 
ready have high-accuracy predictions in the nonlinear regime. 
The idea behind this approach is that a significant amount of 



FIG. 4. — Test of realization scatter on small scales. Shown are the results 
from six simulations that have the same realization on large scales but different 
realizations on small scales. The red line indicates the transition between the 
two regimes. Note that some of the realization noise from the small scales 
leaks back to larger scales. 



timate of the amount of mismatch between the 1300 Mpc and 
365 Mpc Coyote simulations that we can expect after matching 
their power spectra. There is also a small amount of leakage of 
power from small to large scale modes, since the power spectra 
only match at k <C 1 .39 Mpc" 1 , where they were seeded identi- 
cally in the initial conditions. 

To summarize, finite box size effects are clearly an impor- 
tant issue. Some of the problems such as realization scatter 
can be overcome by generating a large number of realizations 
(though this is expensive) but some small inaccuracies in the 
power spectrum will be unavoidable until larger volume, high- 
particle loading simulations can be carried out. 

4. HUBBLE PARAMETER EXTENSION 

In our original work, the value for the Hubble parameter was 
automatically determined from the other five cosmological pa- 



k[1/Mpc] 



FIG. 5. — Accuracy of the linear power spectrum emulator of the full six- 
dimensional parameter space. Shown is the ratio of the emulator to the linear 
power spectra for ten extra models not used for constructing the emulator itself. 
The error overall is within 1%. Note that we only included information on the 
Hubble parameter h on large scales, with an upper cutoff at k = 0.03 Mpc -1 . 
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FIG. 6. — Effects of the six different parameters within their respec- 
tive prior ranges on the linear matter power spectrum. Shown is the ratio 
with respect to a model where each parameter is fixed at its median value: 
9 = {0.0225,0.1375,0.95,0.7,-1.0,0.755} and the sixth parameter (denoted 
at the top of each panel) is varied as a function of k with &o=l Mpc -1 . Light 
blue colors show results for small values of the sixth parameter, while pink 
colors show results for large values. The change of the power spectrum when 
varying the Hubble parameter h is independent of scale, the reason why in the 
case of the linear power spectrum the addition of h as a free parameter works 
well, even when using information only on very large scales. 



information is readily available in the linear to mildly nonlin- 
ear regime covered by RPT. In addition, the very accurate pre- 
dictions at low k, as provided by the 100 power spectra, helps 
to anchor the power spectra correctly on large scales. The in- 
formation on small scales is then provided by the available 37 
models that have been fully simulated. 

We modify the emulation procedure from Lawrence et al.l 
(2010) to account for the inclusion of both 'long' power spectra 
over the entire k range of interest a nd 'short' powe r spect ra in 
the low k range only. Recall that in [Lawrence et al.l (12010). the 
power spectra are modeled using a basis representation: 



a p 



V(k;z;0)=y <j>i(k;z)Wi(0), 6>e[0,lf 
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where the <fii(k;z) are the basis functions, the W0) are the corre- 
sponding weights, and the 9 represent the cosmological param- 
eters. The basis vectors are constructed from principal compo- 
nents and the weights are modeled as a Gaussian process. We 
again use Gaussian processes to model the weights, but con- 
struct the basis vectors slightly differently to better represent 
the variation in the long and short power spectra. 

Two features inform our selection of basis vectors: the in- 
clusion of RPT power spectra that cover only low k and the 
combining of power spectra from A^-body simulations that use 
different box sizes. First, we compute three principal compo- 
nent basis vectors using the 37 smoothed spectra resulting from 
combining A^-body simulations of different sizes. Three basis 
vectors is sufficient to capture the systematic variation at high k 
and the shorter RPT spectra yield well-behaved weights when 
projected on to these. Second, we remove the effects of the first 
three basis vectors and compute three additional PC basis vec- 
tors on the residuals of the 37 smoothed spectra, but only over 



the k range of the original simulations. This improves the re- 
sulting prediction of spectra over this k range, while avoiding 
extra variation at high k resulting from the combination. Again, 
the residuals for the RPT power spectra project on these basis 
vectors in a well-behaved manner. Finally, we compute three 
additional basis vectors on the residuals from all 137 simula- 
tions over the k range covered by the RPT power spectra. This 
further improves the modeling over this low k region. For both 
sets of basis vectors that cover only part of the k range, the basis 
vectors taper linearly to zero over their last 50 k values to en- 
sure continuity. This gives a total of nine basis vectors whose 
weights are modeled with a Gaussian process. 

In order to test the build-strategy that allows us to combine a 
different number of simulations for small and large k, we carry 
out a test with the linear power spectrum. We generate 100 pre- 
dictions for the linear power spectrum out to k < 0.03 Mpc" 1 
and in addition 37 predictions following the original design 
over the full k range out to k < 10.0 Mpc" 1 . We then build 
an emulator over the full k range allowing h to vary. In order 
to test the accuracy of the new emulator, we generate a set of 
10 additional power spectra over the full parameter range and 
compare the emulator prediction with those exact results. Fig- 
ure|5]shows the ratio of the emulator prediction with the linear 
theory answer. Overall, the accuracy is around 1%, which is 
very satisfactory. 

With the new emulator in hand, we briefly study the depen- 
dence of the power spectrum on the different cosmological pa- 
rameters. For this we generate a set of sensitivity plots, shown 
in Fig. [6] The idea behind the sensitivity study is simple: we de- 
termine the power spectrum for a model in the center of the pa- 
rameter hypercube and then vary one parameter at a time from 
its smallest to highest value. This illustrates concisely the ef- 
fect of each parameter on the power spectrum on all scales stud- 
ied. The Hubble parameter primarily shifts the amplitude of the 
power spectrum up and down (right panel, middle row), which 
explains why the simple addition scheme here works so well. 

We also produce a second emulator that keeps h fixed. This 
emulator does not use the additional RPT spectra, but does in- 
clude the results from the smaller A^-body simulations. This 
emulator was constructed in the usual manner, with 6 PC basis 
functions computed from the 37 spectra over the whole k range. 

We will show in a future publication that the idea of com- 
bining different numbers of models at different length scales is 
also very useful in obtaining accurate predictions of the power 
spectrum on even smaller scales. As alluded to in the Introduc- 
tion, computing power spectra on such scales is computation- 
ally very expensive. It is very helpful if the number of models 
needed for this can be kept to a minimum. 

5. SMOOTH POWER SPECTRUM GENERATION 

The power spectra from the A^-body simulations are smoothed 
with essentially the sa me process convolution model used in 
lLawrence et al.l (1201 Ol) . but with an important addition to ac- 
count for vertical shifts where power spectra from different sim- 
ul ati ons are pasted to g ether. 

In lLawrence et al.l (1201 Ol) . the simulation of cosmology c, 
resolution s, and replicate i produces a spectrum, P c si . We 
model this as a multivariate Gaussian variable with a known 
covarianc e 17, and a sm ooth mean described by a process con- 
volution (Hig donll2002l) . A process convolution is constructed 
by generating a latent stochastic process and smoothing it. In 
lLawrence et al.l d2010l) . the spectra for each cosmology share a 
latent process, u c , modeled as a Brownian motion observed on a 
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sparse grid. These latent processes are smoothed by a common 
smoothing matrix, K" ', made from Gaussian smoothing kernels 
whose kernel width varies across the domain in order to ac- 
count for nonstationarity. The resulting model has a probability 
density function 

,/|l/2 
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where the matrix A s truncates the length of the spectrum de- 
pending on the resolu tion. The modeling details are given in 
lLawrence et al.1 (l2010h . 

For the current work, we make a small addition to the mean 
function for the spectrum for each simulation. The spectra for 
the highest resolution simulations are augmented with spectra 
from the smaller boxes as described in the assembly tables. At 
each matching point, the spectra from each box may be ver- 
tically offset because of shot noise. However, we know that 
each simulation is approximating a smooth spectrum across the 
k range. To account for this, we modify the mean structure for 
the simulation model. The mean structure is still built around 
a process convolution model, but includes additional terms that 
move each section of the simulated spectrum up or down. 

Let H c be a matrix with a row for each k value and three 
columns (one for each of the small boxes - no offset if estimated 
between the perturbation theory and the original set of simula- 
tions). Let the cosmology index c also index redshift (although 
unstated, this is also true in Lawrence et al. 2010). This matrix 
describes the matching. At a particular k, if the simulated spec- 
trum is represented by a result from perturbation theory or the 
original simulations, the matrix H c has a row of zeros. If the 
simulated spectrum uses the result from the 365 Mpc box, the 
row has a one in the first column. If the simulated spectrum 
uses a result from the 180 Mpc box or the 90 Mpc box, the row 
has a one in the second or third column respectively. Each can 
have at most a single non-zero entry. Let b c be a set of three co- 
efficients that represent the offset of the simulated result from 
the unknown smooth mean. We add these terms to the model. 
Let V c =K a u c +H c b c , which gives the density 

t'l 1 / 2 
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This model conveys the information that there is a smooth spec- 
trum represented by K a u c and that we observe a noisy version 
of it that has vertical shifts over some ranges. 

The offset parameters need a prior distribution to complete 
the specification. We give them a zero-mean Gaussian prior: 

TT(b)ocexp\~b'b\, (14) 



The parameter A^ is a precision parameter that we set equal 
to the known precision of the simulated spectrum at the first 
matching point. Like the latent process u c , the b can be inte- 
grated out of the posterior and need not be drawn as part of the 
Markov chain Monte Carlo approach. 

6. EMULATOR FOR THE MATTER POWER SPECTRUM 

The construction of the new power spectru m emulator now 
follow s the general methodology laid out in iHeitmann et al.l 
(2009J) and lLawrence et al.l (1201 Oh with the small modifica- 
tions discussed above to combine power spectra over different k 
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FIG. 7. — Emulator predictions for the nonlinear power spectra for the best 
fit cosmologies for WMAP-7 and Planck (see text for details on the parameter 
choices) at z = 0. Due to the higher values for cj„, and a% , the Planck results 
imply a noticeable change in the nonlinear regime. 



ranges. In this section we present some results from the emula- 
tor and tests of its accuracy. To start, we show the predictions of 
the emulator for the best-fit cosmology found by the WMAP-7 
and Planck surveys in Fig. [7] Following the notation in Eq. (01 
the parameters used are 

^wmap-7 = {0.0225,0.13328,0.97,0.7,-1.0,0.81}, (15) 

tfpianck = {0.02225,0.14026,0.968,0.6816,-1.0,0.8284}, 

as obtained from CMB measurements alone. While the spectra 
are close, the Planck parameters lead to an enhancement of the 
power spectrum in the nonlinear regime. 

Below we describe a number of tests of the emulator ac- 
curacy. We find that the emulator quality is better than 5% 
over a large k and z range, even when including the new free 
parameter, h. We also include a comparison against the lat- 
est improved version of Halofit as implemented in CAMB (see 
lTakahashietaill2012l for details). 

6.1. Power Spectrum Predictions 

As explained in Section [2] the emulator is built using results 
from 37 cosmological models, specified in Table Q] In addition 
to these models, we have generated P(k) results for a ACDM 
model (M000), which is close to the current best-fit measure- 
ments. Simulations for this model were carried out in exactly 
the same way as for the other models, using several realiza- 
tions for each box size, and then matching them to produce a 
smoothed power spectrum. This reference P(k) can be used as 
one test of the emulator's accuracy. 

Figure [8] shows results for two emulator versions, one in 
which the value of h is locked by the CMB distance to last 
scatt ering constraint for e ach model, as in the original emula- 
tor (lLawrence et alj|2010b (upper panel), and a second version 
where h is allowed to be a free parameter following the ap- 
proach presented in Section|4](lower panel). Both panels show 
the ratio of the emulator to the simulation result as a function of 
k, the different solid curves correspond to eleven different red- 
shifts between z = 4 and z = used to build the emulator. In the 
/i-fixed case, over a wide k range, and for all redshifts, the em- 
ulator prediction is accurate at the 1 % level, degrading slightly 
only at larger k values, but still remaining better than 5%. 
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FIG. 8. — Accuracy of the emulator predictions for the reference ACDM 
model, M000, not used to construct the emulator. The ratio of the emulator 
prediction to the smoothed M000 simulation result is shown as a function of k. 
In the upper panel we show the predictions for the emulator with the Hubble 
parameter h fixed to the best-fit CMB value for each model, in the lower panel, 
the results when /; is allowed to be an independent parameter. The different 
colors represent results at different redshifts: blue: a = 0.2, green: a - 0.25, 
red: a = 0.2857, cyan: a = 0.3333, purple: a = 0.4, yellow: a = 0.5, gray: 
a = 0.6, blue: a = 0.7, green: a = 0.8, red: a = 0.9, cyan: a = 1.0. The sequence 
of a values correspond to the redshifts, z = 4, 3, 2.5, 2, 1.5, 1, 0.67, 0.43, 
0.25, 0.11, 0.0. The dashed horizontal line indicates the 1% accuracy limit. 
Allowing for h to vary freely only mildly affects the emulator error behavior. 



FIG. 9. — Holdout test for six models. For every test, the chosen model is 
excluded while building the emulator. The emulator prediction is then com- 
pared with that of the 'held out' model. The ratio of the emulator prediction 
to the smoothed simulation result for the six reference models is shown. As in 
Fig. [8] the top panel shows results with ft-fixed emulators, while in the bottom 
panel, this constraint is relaxed. The different colors show results for different 
models (for each model we show results at all redshifts): blue: M004, green: 
M008, red: M013, cyan: M016, purple: M020, yellow: M026. The dashed 
horizontal line indicates the 1% accuracy limit, while the dashed vertical line 
shows k = 1 Mpc~' , the limit of our previous emulator (Lawrence et al. 2010). 



When h is allowed to range freely (lower panel), then over 
the k range where RPT is used (up to k ~ 0.03 Mpc" 1 for the 
low redshift results), the result is excellent, better than that in 
the top panel. This is not surprising since for this range we 
now have 100 models to predict the power spectrum. Beyond 
that matching point, the accuracy degrades slightly in the quasi- 
linear regime compared to the more restricted emulator, but is 
still around 1%. Finally, beyond k ~ 1 Mpc" 1 , the predictions 
are less accurate, but the worst case error remains around 5% . 

Next, we present results from a number of 'holdout' tests. In 
these tests, one model out of the 37 is excluded and an emula- 
tor is constructed based on the remaining 36 models. Then with 
this new emulator, a prediction for the excluded model is gen- 
erated and the accuracy of the prediction determined. This test 
has an obvious caveat, particularly relevant in cases where not 
many models are available - degradation of the emulator qual- 
ity because an important point in the parameter-space hyper- 
cube is omitted. This can be particularly serious if the excluded 



point is on the edge of the design, because now extrapolation to 
the edge of the hypercube is required, which has a much higher 
error. In order to avoid this additional complication, we restrict 
our holdout test to models that lie well within the hypercube. 
Nevertheless, it should be kept in mind that reported errors in 
the holdout tests can be larger than the actual values. The hold- 
out test results, shown in Fig.|9]are similar to the ones found for 
M000, both in error amplitude and trends with increasing k. In 
the case where h is allowed to be a free variable, the predictions 
are again less accurate, already in the quasi-linear regime, and 
degrade further beyond k ~ 1 Mpc" 1 , to of order 5%. Note that 
the low-A: error band is consistent with that for the linear theory 
test case (Cf. Fig . [5j ■ 

To summarize, the accuracy of the new, extended, emulator 
is well-characterized at all redshifts. With h fixed to the CMB 
constraint, the accuracy to k ~ 1 Mpc" 1 is at the ~ 1% level, 
degrading to ~ 5% for k ~ 10 Mpc" 1 . When h is allowed to 
be a free parameter, the degradation in accuracy is relatively 
modest. 
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FIG. 10. — Same as in Fig. [6] but for the nonlinear power spectrum. The 
upper panel shows results at z = 4 while the lower panel shows results at z = 0. 



6.2. Sensitivity Investigation 

With the full emulator at hand we can now carry out a sensi- 
tivity analysis for the nonlinear power spectrum, along the same 
lines as conducted in Section|4]for the linear power spectrum. 
As before, we fix five of the six cosmological parameters at the 
midpoints of their range and then vary the sixth parameter be- 
tween its minimum and maximum value. The results - for the 
lowest and highest redshift we consider - are shown in Fig. [10] 
The upper six panels show the results at z = 4 and the lower six 
panels at z = 0. 

Some interesting features in the results can be noted. As to be 
expected, the baryon fraction wj does not affect the power spec- 
trum significantly. The best constraints we have for Ub come 
from CMB measurements; Planck determines this parameter to 
exquisite accuracy. The effect of the matter fraction ui m domi- 
nates at large scales, here again, CMB results deliver very good 
estimates of this parameter. In general, it is important to allow 
for variations in u,„ due to degeneracy issues; u m influences 




FIG. 11.- — Performance of the emulator in comparison with Halofit from 
Takahashi et al. 12012) for M000 at z = 0. Shown are the ratios of the emula- 
tor with fixed Hubble parameter and the smoothed simulation (blue), with the 
free Hubble parameter emulator (red), and with Halofit (black). (The emulator 
results are shown also in Fig.[8]in cyan and are here repeated for easy compar- 
ison with Halofit). The A-fixed emulator is accurate at the percent level out to 
k ~ 1 Mpc~' and at the 2% level at higher k, the ft-free emulator is accurate at 
the 2% level throughout. Halofit on the other hand shows deviations at the 3% 
level at k ~ 0. 1 Mpc~' and up to 6% in the higher k regime. 



the overall amplitude of the power spectrum as do as and h. 
The spectral index n s , influencing the tilt of the power spectrum 
also mainly effects large scale behavior. The remaining three 
parameters, h, w, and in particular a%, alter the power spectrum 
on all scales and constraints on these parameters can therefore 
be improved by measurements of nonlinear scales. 

In particular, the influence of a% at the two redshifts shown 
in Fig. [10] on the power spectrum is noteworthy. At high red- 
shift, nonlinearities for models with large a% (pink) have al- 
ready developed considerably, therefore the ratio with the mid- 
point power spectrum shows large values. The power spectra 
for small values of erg (blue) show only mild nonlinear growth 
and therefore the ratio is still almost flat. Over time, the non- 
linear turn-over moves in to smaller k ranges (affecting larger 
and larger scales), causing the bump just before k = 1 Mpc" 1 
shown in the panel for z = 0. At this epoch the overall differ- 
ence between the models has decreased, since at this point all 
models have reached the nonlinear regime. A similar observa- 
tion can be made for w - at early times the differences between 
the models are more pronounced than at later times. It is easy to 
imagine using the new emulator as a convenient tool to investi- 
gate the effects of varying cosmological parameters and redshift 
on P(k) in the nonlinear regime. 

6.3. Comparison with Halofit 

We now perform a comparison between the emulator and 
Halofit, a popular fitting formula for the matter power spec- 
trum motivated by the halo model. We use the updated version 
with 35 fitting parameters provided by Takahas hi et al.l (1201 2|) . 
which is calibrated to larger volu me and higher reso lution N- 
body simulations than the original ISmith et all (12003) formula. 
It is based on simulations run on six WMAP cosmologies (years 
1, 3, 5, 7, and two WMAP7 models, but with w = -0.8 and 
w = -1 .2, instead of w = -1). The obtained fitting formula was 
checked against simulations run using the first 10 of the 37 Coy- 
ote models. This improved version of Halofit is stated to be ac- 
curate to ~ 5% in the range < k < 1 /zMpc" 1 and ~ 10% at 
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FIG. 12. — Performance of the emulator in comparis on with linea r theory 
(dotted) supplied by CAMB and Halofit from Takahas hTeTaT] 42012D (solid). 
Unless otherwise stated in the legend, the cosmological parameters are: Lo m = 
0.1296,aj 6 =0.0224,fi., =0.97,ct 8 = 0.8,w = -1,/i = 0.72. The parameters for 
WMAP-7 and Planck are given in Eq. )15> . See text for discussion. 



1 < k < 10 /iMpc" 1 , for < z < 3. 

One important feature of the emulator-based approach is that, 
close to the parameter values of the underlying simulations, the 
emulator error is small, if not essentially the simulation error. 
However, this is not the case with a fitting formula. For in- 
stance, even though the N-body simulations that underlie the 
new Halofit formula agree with the original Coyote emulator at 
the ~ 1— 2% level over its ra nge of validity for the 10 Coyote 
models used as a reference in Takahas hi et al.l (12012) (se Fig. 2 
of the cited reference), this good level of agreement drops to 
8-13% (as a function of the k range) when the fitting formula is 
used in the comparison. Depending on what one may be trying 
to achieve, this loss of accuracy in Halofit is significant. 

We first compare the predictions of the emulator (both ver- 
sions, h free and h fixed) at z = for M000 against the smooth 
simulations (in the same way as shown in Fig. [8} and the predic- 
tion of Halofit for the same model. Since the emulator was built 
without including this simulation, this is a good test of both pre- 
diction schemes, emulator and Halofit. The results are shown in 
Fig- HH The emulator predictions for the /j-fixed version are ac- 
curate at 1% out to k ~ 1 Mpc" 1 and then degrade very slightly 
to 2%. The /z-free version of the emulator is accurate at the 1- 
2% level throughout. Halofit shows 3% deviations compared to 
the simulation in the mildly nonlinear regime (k ~ 0.1 Mpc" 1 ) 
and more than 6% in the nonlinear regime. 

Next we show a comparison of the emulator (/z-free version 
only) directly with Halofit for a variety of cosmologies at z = 0. 
Figure Q~2] displays the ratio of the emulator with respect to 
linear t heory power spectra fro m CAMB and Halofit as mod- 
ified bylf akahashi et all (120121) . We choose to test a WMAP7 
dKomatsu et alj|201 lh and a Planck dAdeetal]|2013l) cosmol- 
ogy and two other cosmologies on the limits of the design of our 
emulator. The red and blue curves in particular are on the edge 
of the new design with has a free parameter. The cosmological 
parameters for these two model are the same as model M000: 
uj m = 0.1296, w b = 0.0224,n, = 0.97, cr 8 = 0.8, w = -1,A = 0.72 
and in our comparison we vary one to two parameters at a time 
as noted on the legend. 

At large scales, k < 0.01 /iMpc" 1 , all power spectra are con- 
sistent at the percent level; this is less trivial than it seems since 



Halofit and linear theory start diverging at k ~ 0.003 hMpc~ l 
and we use the full version of the emulator with extreme values 
of h located at the edge of the design. On quasi-linear scales, 
there is ~ 5% deviation with some oscillatory structure from 
the BAO feature. On small er scales, the a ccurac y is consis- 
tent with what is stated in iTakahashi et al.l d2012l) . apart from 
the matter dominated case with il m = 0.51, in which the Halofit 
power spectrum is suppressed by more than 20%. The relatively 
uniform error behavior of an emulator - a result of the sampli ng 
theory and the Gaussian process based-regression methodology 
- is hard to reproduce in a fitting formula; the last result of the 
comparison provides an example of this. 

7. CONCLUSION AND OUTLOOK 

High accuracy predictions for cosmological observables in 
the nonlinear regime will be crucial in the future to exploit the 
power of ongoing and upcoming cosmological surveys. These 
predictions have to at least match, better yet, exceed the accu- 
racy of the measurements themselves. Achieving this goal with 
first principles predictions or perturbation theory is impossible 
because of the highly nonlinear and dynamical nature of the 
problem. High-accuracy predictions must therefore be obtained 
from state of the art simulations. 

In a series of recent papers (the Coyote Universe papers) we 
have shown how to build a prediction tool from a relatively 
small set of simulations (37 cosmological models) to generate 
the matter power spectrum out to k ~ 1 Mpc" 1 at the 1% accu- 
racy out to z = 1 • In the current paper, using what we consider 
the 'minimal' amount of work, we have extended the emulator 
to significantly higher values of redshift and wavenumber. In 
addition, we have added one new free parameter, h. 

With current computational resources, it is impossible to ob- 
tain high emulation accuracy out to large k with a collection 
of single-box simulations, because the dynamic range require- 
ments impose a very high computational cost. The use of nested 
boxes avoids the computational cost, at the expense of reduced 
accuracy, resulting from increased sampling variance as well as 
reduced accuracy in the region of the nonlinear turnover. Nev- 
ertheless, with this approach we have improved on the accuracy 
attained by any other existing prediction tool. Overall, our pre- 
dictions for the power spectrum are accurate to better than 5%, 
and over smaller k ranges around 1%, for < z < 4. In the 
high-fe regime, baryonic effects will have to be included in any 
case, so any future strategy should include methods for improv- 
ing the accuracy of the gravity-only simulation results, as well 
as ways to model the baryonic contributions. 

In a companion paper we will present a method for extrapo- 
lating the results beyond k ~ 10 Mpc -1 (including estimates of 
the associated errors) and an emulator for the shear power spec- 
trum and other weak lensing observables. More detailed work 
relevant for surveys (DES, LSST) is also underway. 

For future surveys, the cosmological parameter space will 
have to be enlarged beyond what was considered here. In the 
case of DES, for example, predictions for dynamical dark en- 
ergy models are needed. Increasing the sample space will fur- 
ther increase the associated computational costs. On the pos- 
itive side, results from surveys such as Planck help to narrow 
down the parameter ranges substantially, which in turn will help 
to reduce the number of models we have to investigate. Multi- 
level sampling schemes can be devised to deal with these situa- 
tions, one of our directions for future work. 

APPENDIX 
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In the course of the work presented in this paper, some 
chan ges were made to pre vious versions of the CosmicEmu 
from Lawrence et al. (2010). 

First, Taruva et al. (2012) present work on perturbation the- 
ory and compare their results to the CosmicEmu for the cos- 
mologies that were simulated to build the CosmicEmu. Over- 
all, they note a close match, but for M015, they find a notable 
difference. In the course of investigating this issue, we discov- 
ered a typographical error in the input files that were used to 
build the CosmicEmu. The file containing the input parame- 
ters differed from the input parameters used to run the actual 
simulations for three cosmologies: M015, M017, and M019. 
The effects of this error were quite local to the neighborhood of 
these three cosmologies. 

Second, we experimented with the prior parameters for the 
Gaussian process emulation procedure. In order to explain this 
more fully, we must describe o ne more subtle d etail in Gaussian 
process fitting. Appendix B in lHeitmann et alj J2009) describes 
the statisti cal model for the Gauss ian process emulation. Equa- 
tion B 1 in iHeitmann et alj (|2009) gives the Gaussian process 
distribution for the principal component weights with covari- 
ance matrix £ = X~ l R(6,9';p Wi ). The function R(-), given in 
Equation B2, results in a response that interpolates and is ex- 
tremely smooth. As a result, the model can be susceptible to 
estimation problems and overfitting. Thus, in practice, we in- 
clude an additional term in the covariance specification to relax 
the interpolation requirement slightly. The actual covariance is 
given by E = X^ t R(9,6';p Wl ) + X^l g I, where I is the identity ma- 
trix. The new parameter \ mlg is an additional precision term 
(the 'nugget') that governs how well the estimated response 
function interpolates the data. When this parameter is large, the 
response surface interpolates better. We estimate this parame- 
ter along with all of the others, so the data can inform about the 
smoothness of the model. Like the other precision parameters 
described in Heitmann et al. (2009), it has a Gamma prior dis- 
tribution (see Equation B5 for an example). Unfortunately, the 
default prior parameters in the estimation software can prevent 
these precisions from being as large as the data would allow and 
prevent the response surface from capturing all of the behavior 



of the simulations. We now set these parameters at a 
and b 



nug,i 



nug,i ' 



1 

where i indexes the principal components. This 
prior gives the data much more control over the parameter es- 
timate. The new estimates for X nug are significantly larger than 
the previous estimates. The resulting response surface for each 
principal component basis is more accurate and more principal 
component basis functions can be used (seven instead of five). 
Figure[T3]shows the results of these two improvements on the 
holdout predictions and the prediction for the best fit cosmology 
M000. The fits, while good before, are now further improved. 
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FIG. 13. — Holdout (top) and M000 (bottom) results for th e improved fit 
to the original CosmicEmu described in Lawrence et al. (2010). Compare to 
Figures 9 and 10 in that paper. 
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